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> ABSTRACT 

o 

^7' Context. Using the recent observational data of Roser et al. we present TV-body simulations of the Hyades open cluster. 

Aims. We make an attempt to determine initial conditions of the Hyades cluster at the time of its formation in order 
CO to reproduce the present-day cumulative mass profile, stellar mass and luminosity function (LF). 

Methods. We performed direct TV-body simulations of the Hyades in an analytic Milky Way potential that account for 
' , ' stellar evolution and include primordial binaries in a few models. Furthermore, we applied a Kroupa (2001) IMF and 

used extensive ensemble-averaging. 

Results. We find that evolved single-star King initial models with King parameters Wo = 6 — 9 and initial particle 
numbers A^o ~ 3000 provide good fits to the observational present-day cumulative mass profile within the Jacobi 
radius. The best-fit King model has an initial mass of 1721 Mq and an average mass loss rate of —2.2 TVf0/Myr. The 
^ I K-band LFs of models and observations show a reasonable agreement. Mass segregation is detected in both observations 

and models. If 33% primordial binaries are included the initial particle number is reduced by 5% as compared to the 
model without primordial binaries. 

Conclusions. The present-day properties of the Hyades can be well reproduced by a standard King or Plummer initial 
model when choosing appropriate initial conditions. The degeneracy of good-fitting models can be quite high due to 
the large dimension of the parameter space. More simulations with different Roche-lobe filling factors and primordial 
binary fractions are required to explore this degeneracy in more detail. 
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1. Introduction Portegies Zwart et al. (2001) differs from the present study 

by the fact that they did not include a kick for the white 

The Hyades (the "rain stars") are one of the best-known dwarfs but argued that they could be hidden in binary sys- 

open star clusters on the northern sky. They are located in tems with considerably brighter companion stars. Madsen 

the constellation Taurus,^ close (on the celestial sphere) to (2003) concentrated on investigating the accuracy of as- 

thc Pleiades or "Seven Sisters" , the most prominent open trometric radial velocities and compared the HR diagrams 

star cluster on the northern night sky. The Hyades and the and internal velocity dipersions of observations and models. 

Pleiades form the "golden gate of the ecliptic" , i.e. the eclip- On the other hand, Chumak, Rastorguev & Aarseth (2005) 

tic separates the two clusters on the celestial sphere. The studied the Hyades orbit and provided first models of the 

C3 Hyades cluster is markedly older than the Pleiades clus- tidal tails of the Hyades. 

ter. It has an age of 625 ± 50 Myr derived from the helium xhe starting point for the present work is an obser- 

abundance in combination with isochrone modeling which yational data file (Roser et al. 2011) from which masses, 

includes convective overshooting (Perryman et al. 1998). positions and velocities for 724 probable member stars 

The shape of the Hyades Hertzsprung-Russell (HR) dia- of the Hyades can be derived. The PPMXL catalogue 

gram with a narrow and well-defined main sequence has (Roser, Demleitner & Schilbach 2010) contains positions 

been recently studied by von Leeuwen (2009). Previous and proper motions of the Hyades members down to 

A^-body simulations of the Hyades have been performed 0.II6 MqQ The membership has been determined with 

by Portegies Zwart et al. (2001), Madsen (2003; see also the convergent point method: Given the right ascension, 

references therein) and Chumak, Rastorguev & Aarseth declination and the respective proper motions, the conver- 

(2005). These studies are in many respects similar to the 

present study (although special emphasis is placed on dif- ^ Note that there is a recent study (Goldman et al. 2011) based 

ferent aspects of the modeling of the Hyades). The study by on PanSTARRSl which goes down to even lower masses. 
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Table 1. Hyades parameters as determined from the ob- 
servational data file. For the determination of the Jacobi 
radius we assumed that the solar radius is at Rq = 8 kpc 
(values for Rq = 8.5 kpc are given in brackets). 



Quantity 


Value 


Lowest mass rrii 


0.116 Mq 


Highest mass 


2.6 Mq 


Maximum radius rmax 


30 pc 


Number of stars A'^(30 pc) 


724 


Total mass M (30 pc) 


469 Mq 


Jacobi radius rj 


8.6 pc (9.0 pc) 


Number of stars N{rj) 


354 (359) 


Enclosed mass M{rj) 


272 M0(275 Mq) 


Core radius Vc 


3 pc 


Number of stars A*' (3 pc) 


97 


Enclosed mass M(3 pc) 


99 Mq 


Completeness limit 


0.25... 0.17 Mq 





D o 1, Total 

Hf, = a Kpc g^ilqe 




Vo = 234.2 km/s Disk 




Halo 















Fig. 1. Rotation curve (at z = 0) of the three-component 
Plummer-Kuzmin model of the Milky Way. 



Table 2. The list of galaxy component parameters. The 
first column gives the component, the second the mass, 
and the third and fourth the Plummcr-Kuzmin parameters 
(equation A.l). 



Component 


M [Mq] 


a [kpc] 


b [kpc] 


Bulge 


1.4 X 10"' 


0.0 


0.3 


Disk 


9.0 X 10^" 


3.3 


0.3 


Halo 


7.0 X 10" 


0.0 


25.0 



X [kpo] 

Fig. 2. "Artistic" RGBA image of the (logarithmic) density 
structure of the three-component Plummer-Kuzmin model. 
The three components are overblended with an alpha chan- 
nel with 90% transparency. The Halo is plotted in blue, the 
Disk in green, and the Bulge in red. 



gent point method predicts for each candidate member a 
heliocentric distance (secular parallax) and radial velocity. 
Thus the full 6D phase space information is given. Roser 
et al. determined the individual masses from the mass- 
to-luminosity relation by Pinsonneault et al. (2004) and 
from Dartmouth isochrones (Dotter et al. 1998) and BCAH 
isochroncs (Baraffe et al. 1998). Binarity is known only for 
a minor portion of the sample (Roser et al. 2011). 

Table [T] gives a few parameters of the Hyades cluster as 
determined from the observational data file. We determined 
the tidal (or Jacobi) radius (King 1962) iteratively from the 
observational data according to the procedure described in 
Section 4 of Ernst et al. (2010) assuming that the solar 
radius is at _Ro = 8 kpc (the values for Rq = 8.5 kpc are 
given in brackets) and that /3 = n/fl = 1.37, where k and 
are the epicyclic and circular frequencies of a near-circular 
orbit at the solar radius. The incompleteness in the stellar 
mass function sets on somewhere between 0.25 and 0.17 M0 
(Roser et al. 2011). 

This paper is organized as follows: Section [2] describes 
the numerical methods which we use. Section (S] describes 
our integration of the Hyades orbit. In Section |4] we discuss 
in detail our parameter space. Section [5] summarizes the 
results. Section [6] contains the discussion and conclusions. 

2. Numerical method 

We solve the A^-body problem of the Hyades in an ana- 
lytic background potential of the Milky Way. For the back- 
ground Milky Way potential we use an axisymmetric three- 



component model, where the bulge, disc, and halo are de- 
scribed by Plummer-Kuzmin models (Miyamoto & Nagai 
1975) with the potential 



y^i?2 + (a + + z2)2 

The parameters a, b, and M of the Milky Way model are 
given in Table [2] for the three components. In the limit 
a — this model reduces to a Plummer model (Plummer 
1911). On the other hand, in the limit 6-^0 the model 
reduces to a Kuzmin model (Kuzmin 1956). 

Fig. [1] shows the rotation curve of the three-component 
model of the Milky Way. As in our previous works, the 
parameters of the three-component model are chosen such 
that the rotation curve matches that of the Milky Way 
(Dauphole & Colin 1995). At the solar radius Rq = 8.0 
kpc, which we have assumed in this study, the value of the 
circular velocity is Vq = 234.2 km/s. The values of Oort's 
constants A and B are consistent with the observed values 
{A,B) = (14.5 ± 0.8,-13.0 ± 1.1) km/s/kpc derived by 
Piskunov et al. (2006). 

Figure [2] shows an "artistic" RGBA image of the (log- 
arithmic) density structure of the three-component 
Plummer-Kuzmin el. The three components are 
overblended with an alpha channel with 90% trans- 
parency. The Halo is plotted in blue, the Disk in green, 
and the Bulge in red. 

For the solution of the A^-body problem in the ana- 
lytic background potential of the Milky Way, the new di- 



2 



A. Ernst et al.: Simulations of the Hyades 



rect A^-body code nbody6tidgpu (see Appendix A for a 
description) is used. Since this program is based on the 
code NBODY6 by Aarseth (see Aarseth 1999, 2003), it is 
possible to follow stellar evolution within the Hertzsprung- 
Russell diagram. The implementation of stellar evolution 
for individual stars of the A^-body system is based on an- 
alytic formulae from which radius, luminosity and stellar 
type as a function of the initial mass and age can be de- 
rived (Aarseth 2003; for the stellar evolution recipes see 
Hurley, Pols & Tout 2000 and Hurley, Tout, Aarseth & 
Pols 2001). In particular, this approach includes the forma- 
tion of stellar remnants like white dwarfs (WDs), neutron 
stars (NSs) and black holes (BHs). In addition to the mod- 
elling of standard evolutionary mass loss by stellar winds 
(which is implemented in many TV-body codes available to- 
day), Aarseth's code and its variants include routines for 
the treatment of large instantaneous mass loss due to spe- 
cial events like the occurence of supernovae and the for- 
mation of planetary nebulae. (We remark that it is not so 
easy to correct energy, all forces and first derivatives within 
the A'^-body system for the case of large instantaneous mass 
loss of one particle in the case of, e.g., a supernova.) Thus it 
is possible with nbody6 (and its variants) to model a star 
cluster realistically even in the first 55 Myrs of evolution in 
which the supernova events occur (e.g. Hansen & Kawaler 
1994, their Equation (1.88)). 

The code nbody6 and its variants optionally apply ve- 
locity kicks to stellar remnants when a supernova occurs or 
a planetary nebula forms. The random kick velocities for 
supernova kicks are drawn from a Maxwellian distribution 
with a ID velocity dispersion of ~ 190 km/s correspond- 
ing to a mean of « 300 km/s (Hansen & Phinney 1997). 
For all WDs, the ID dispersion of the Maxwellian is taken 
to be 5 km/s corresponding to a mean of w 8 km/s (e.g. 
Fellhauer et al. 2003 for a lower limit on that dispersion). 
We remark that the final number of WDs in our models 
depends critically on this dispersion. The escape velocity 
from the center of the cluster is « 6 — 7 km/s for an open 
cluster with a mass of 2000 M© and a half-mass radius of 
3.5 pc assuming that it has a Plummer profile (Plummer 
1911). 

Finally, we note that in nbody6tidgpu, all stars are 
kept in the simulation forever, i.e. there is no optional re- 
moval of escapers as in nbody6 and its other variants. The 
reason is that we would like to follow the stellar orbits in 
the tidal tails and investigate their properties. 



3. Orbit 

We have integrated the center of mass orbit of the Hyades 
backwards in time in the analytic Milky Way potential (i.e. 
the three-component Plummer-Kuzmin model described in 
section [2]) . We neglected encounters with giant molecular 
clouds and the effect of spiral arm passages and disk shock- 
ing. The integration time (625 Myr) was equal to the most 
probable age of the Hyades (Ferryman et al. 1998). The full 
3D orbit is shown in Figure [3] The orbit integration with a 
simple integrator yielded the initial position of the Hyades 
at the time of its formation. From the initial position of the 
Hyades we ran full A^-body models up to the present time. 

The present-day position of the Hyades was taken to be 



xo = -44.26 pc - 8000.00 pc, (2) 
Va = 0.31 pc, (3) 
zo = -16.89 pc + 20.00 pc. (4) 

The first terms in Equ. ([2]) - Q correspond to the center 
of mass of the observed Hyades spatial distribution with 
respect to the position of the Sun. The second terms in 
Equ. ([2]) and Q represent the solar position in the Milky 
Way. We used (i?o, -zo) = (8.00,0.02) kpc (consistent with 
Piskunov et al. 2006). 

The present-day velocity was given by 

vm = 41.18 km s"^ - 11.10 km s"\ (5) 
VyQ = 19.04 km s"^ - 12.24 km s"^ + 234.20 km s"\ (6) 
w^o = 1-27 km s"^ - 7.25 km s"\ (7) 

The first terms in Equ. ([s]) - ([T]) represent the velocity of 
the center of mass of the observed Hyades distribution with 
respect to the Sun. The second terms in Equ. ([5]) - ([t]) 
correspond to the peculiar motion {U,V,W)q of the Sun 
with respect to the local standard of rest (LSR) according 
to Schonrich et al. (2010). The third term in Equ. (pi) is 
the orbital speed of the LSR at i?o = 8 kpc in the three- 
component Plummer-Kuzmin model of Section [2j 

Note that with these parameters the Hyades orbit 
around the Galactic center is oriented in a clockwise sense 
as seen from the Galactic north pole. The orbit has a con- 
siderable vertical oscillation in 2;-direction between ss ±74 
pc. In the radial direction the orbit is mildly eccentric and 
oscillates between 7087 pc and 8638 kpc. Roughly three 
orbits are completed in 625 Myrs. 

4. Parameter space 

We were looking for the best fitting King models (King 
1966) and Plummer models (Plummer 1911) to the obser- 
vations of Roser et al. after 625 Myrs of evolution with 
respect to the final cumulative mass profile M[r) within 




Fig. 3. 3D orbit of the Hyades in the three-component 
Plummer-Kuzmin model of the Milky Way. The integra- 
tion time is 625 Myrs. 
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< r < 9 pc (Jacobi radius). We compare the cumula- 
tive mass profiles of models and the observational data and 
not the density profiles because the cumulative mass pro- 
files are of a better quality with respect to noise. Secondly, 
we note that the observational data of Roser et al. contains 
only 724 stars but the observational present-day mass func- 
tion (PDMF) drops below the completeness limit which lies 
between 0.25 and 0.17 Mq. If we assume that the PDMF 
is shallow at the low mass end according to Equation [9] 
below, ~ 100 stars would be missing in the low-mass bins 
above 0.116 Mq (see Roser ct al. 2011, their Figure 10). We 
took that into account for the comparison of the particle 
numbers and total masses of observation and model. 

All our models are initially Roche-lobe filling, i.e. the 
initial 99% Lagrange radius r9g% was set to be equal to 
the tidal (or Jacobi) radius rj. The 99% Lagrange radius 
is defined as the radius which contains 99% of the cluster 
mass. The Jacobi radius is given by (King 1962) 



GM 



-,1/3 



(8) 



where G, M, f2 and k are the gravitational constant, the 
cluster mass, the circular and the epicyclic frequencies of a 
near-circular orbit, respectively. 

We define the initial mass function (IMF) as the number 
of stars per unit logarithmic mass interval. It is given by 



dN 
din m 



ci 

C2 m" 
C3 rrf 
C4 m' 



-0.7 
-0.3 



0.01 < m/Af© < 0.08 
0.08 < m/Af© < 0.50 
0.50 < tti/Mq < 1.00 
1.00 < m/A/0 



(9) 



where ci — C4 are four constants (which are determined by 
the normalization and three conditions for the transitions 
between the mass regimes) . The exponents in our equation 
(|9| are the exponents ao — 03 given in Kroupa (2001), his 
equation (2), adjusted for logarithmic mass bins. Thus our 
IMF is a Kroupa (2001) IMF. We force the lowest and the 
highest mass to take on the values mi = 0.08 Mq and ruh = 
100.00 Mq, respectively. Thus the first power law in ^ 
for the brown dwarf regime is obsolete. In NBODyBtidgpu, 
the individual masses are randomly drawn from the above 
distribution ^ using a routine by Kroupa with a correction 
by Weidner. 

For all models we switched on stellar evolution in 
nbody6tidgpu. For the metalhcity we used the observed 
value for the Hyades Z 0.024 (Ferryman et al. 1998). 

More than 400 individual models have been computed. 
The run time for one model on a dual-Xeon 3.2 GHz with a 
GeForce 8800 GTS 512 graphics card was about 45 minutes. 

In a first step, a grid of 30 King models (King 1966) 
in the parameter space (Wq, N), where Wq is the King 
parameter and N the initial particle number, has been 
integrated in order to get an overview of the evolution 
of different models. We used for the King parameters 
Wq = 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 and particle numbers N = 
3000, 4000, 5000. Thus the grid ceU size in the parameter 
space is given by {AWo,AN) = (1,1000). We did not in- 
clude primordial binaries in these models. 

Using this first grid, we were able to notice trends with 
respect to the choices of N and Wq- The variation in Wq 
turned out to be weak. Also, the grid was too coarse to 
reach the desired accuracy in N. 



In the second step, we therefore continued with an it- 
erative search in the parameter space {Wo,N) for the best 
fitting models to the observations without primordial bina- 
ries. For this purpose we refined our first grid in the pa- 
rameter space with respect to iV to a grid cell length of 
AiV — 125 particles. On the other hand, it turned out to 
be purposeful to coarsen the grid with respect to Wq to a 
grid cell height of AWq = 3, i.e. we restricted the investiga- 
tion to the models with Wq = 3, 6, 9, 12. Thus the grid cell 
size of the second grid in the parameter space was given by 
(AWo,AiV) = (3,125). 

Our criterion for the best fit is that the inner cumulative 
mass profiles within r = 9 pc (Jacobi radius) of observations 
(Roser et al. 2011, their Figure 6) and models (ensemble 
mean, see below) show the best agreement. 

We did not run all models in this second grid. Rather, 
we started from the best fitting model in the first grid for 
a given Wq and varied the particle number N in steps of 
AA''=125 particles in order to find the best fitting model in 
the second grid. 

At the same time, we performed for each initial model 
in the second grid n = 15 runs with different random num- 
ber seeds. From these 15 runs we obtained physical quanti- 
ties by ensemble averaging. For example, we compared the 
observations with the mean cumulative mass profile of 15 
models in order to find the best fit acording to the criterion 
mentioned above. The random number generator is invoked 
mainly for the following purposes: 

— Drawing of initial masses, positions and velocities from 
the given probability distributions (Kroupa 2001, King 
1966), 

— Assignment of kick velocities for WDs, NSs and BHs. 

— Assignment of perihelion, node and inclination of pri- 
mordial binaries. 

We remark that we did not simply ensemble-average the 
initial model from 15 initial models which differed by their 
random number seeds. In particular, we did not perform 
only one single run from such an ensemble-averaged initial 
model but really performed 15 individual runs because 

— we are looking for standard deviations in the physical 
quantities, 

— the modelling of the population of stellar remnants 
should be as realistic as possible. 

In order to verify that we found indeed the best fitting 
ensembles, we looked at the cumulative mass profile of at 
least one of the neighboring ensembles with respect to N. 

The ensembles in the second parameter grid are named, 
for example, en3000W6 {Nq = 3000,1^0 = 6), i.e. the let- 
ters "en" for "ensemble" are followed by the initial particle 
number Nq and the King parameter Wq. 

We also looked for the best fitting ensemble of single- 
star Plummcr models (Plummer 1911) with respect to the 
inner cumulative mass profile (within r = 9 pc) using also 
steps of AA^ — 125. Note that the choice of rg()%/rj fixes 
the Plummer radius while the King models still have the 
free parameter Wq. The Plummer model runs are discussed 
m Section [5231 

We ran the best fitting ensemble of single-star Plummer 
models with primordial binaries in order to see their effect 
on the evolution. For models with primordial binaries we 
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defined the primordial binary fraction and the fraction 
of stars in binaries fsb as 



fsb — 



2Nb 



(10) 



where Ns and Nb are the initial numbers of binaries and sin- 
gle stars, respectively (e.g. Kroupa 1995, JahreiB & Wielen 
2000, Heggie, Trenti & Hut 2006). The period-generating 
function for the primordial binaries was given by (Kroupa 
1995, Eq. lib) 



logio^(^) 



logio Pmin + [expi2X/rj) - 1]}^/" (11) 



with a uniform random variate X G [0, 1]. The parameters 
are 77 — 2.5, S = 45.0, and Pmin = 5 days. For the binary 
eccentricities a thermal distribution has been adopted with 
the eccentricity-generating function = X where e is the 
eccentricity and X G [0, 1] a uniform random variate. We 
did not assume any correlation between the masses of the 
binary components (e.g. Eggleton, Fitchett & Tout 1989). 

5. Results 

For the statistical analysis, we have applied the quantities 



i=l 



1 ^ 

1=1 



(12) 

where n = 15 is the number of models in an ensemble, Xi 
are the individual measurements, x is their mean, a is the 
standard deviation (i.e. the mean error of a measurement 
in a single model) and a„i is the standard error (i.e., the 
mean error of the mean). 

5.1. Cumulative mass profiles 

The contamination with field stars in the observational data 
file is negligible for r < 9 pc, 7.5% at 9 pc < r < 18 pc, 
and 30% at 18 pc < r < 30 pc (with respect to the number 
of stars; Roser et al. 2011), where r is measured from the 
center of the cluster. 

Three comparisons of the cumulative mass profiles are 
shown in Figure |4] The top panel shows the comparison for 
the ensemble of models en3000W6 {Nq = 3000,1^0 = 6). 
Shown is the mean cumulative mass profile of the ensemble 
of 15 models (thick solid red line), i.e. we averaged over the 
15 single cumulative mass profiles of the individual mod- 
els. Compared with the observations (dashed blue line), the 
ensemble en3000W6 is the best fitting of all our models 
according to the criterion stated in Section |4] Assuming 
that the model is perfect, the field star contamination is 
|M^odei-Mobs|/Mobs = 2% (3 pc), 2% (9 pc), 12% (18 pc), 
21% (30 pc), where M^^d and Afobs are theoretical and ob- 
served cumulative masses. 

The second-best fitting ensemble of King models 
(en3000W9; iVo = 3000, Wq = 9) is shown in the middle 
panel of Figure [4j The ensemble mean of the cumulative 
mass profile fits the inner part of the observational pro- 
file also fairly well. The field star contamination is here 
|M^odei-Mob.|/Mobs = 12% (3 pc), 5% (9 pc), 18% (18 pc), 
26% (30 pc). 



Table 3. Parameters of the best fitting ensemble of models 
en3000W6. Notes: (1) The final parameters are obtained 
after 625 Myrs of evolution, i.e. at the present time. (2) 
The final "observed" particle number is the number of stars 
with m > 0.116 Mq. (3) The highest mass is obtained by 
excluding BHs. 



if 


Quantity 




Value 


1 


Initial model (initial) 








Particle number A^o 




3000 




Total mass (Mo) 




1721 M© 




King parameter Wo 




6 




Jacobi radius rj 




16.2 pc 


2 


r < 30 pc (final) 








Particle number {N{} ± o"™ 




736 ± 50 




"Observed" part. numb. (Af)obs 




625 ± 42 




Lowest mass mi 




0.08 Mq 




Highest mass ± ct™ 




2.68 ± 0.01 M0 




Mean mass (m) 




0.503 M0 




Total mass M; ± a^n 




369 ± 25 Mq 


3 


r < 9 pc (final, inside Jacobi radius) 






Particle number (Nc) ± am 




476 ± 43 




"Observed" part. numb. (Af)obs 


± am 


414 ± 37 




Mean mass (m) 




0.565 Mq 




Total mass Mf ± am 




269 ± 25 Mq 


4 


r < 3 pc (final, inside core) 








Particle number (Af) ± am 




124 ± 14 




"Observed" part. numb. (Af)obs 


±am 


112 ± 12 




Mean mass (m) 




0.784 M0 




Total mass Mf ± am 




97± IIM0 



Table 4. 3D velocity dispersions at t = 625 Myr for the 
ensemble en3000W6 for stars within a radius of 30, 9 and 
3 pc from the cluster center. 



Radius 


(aso) ± am [km/s] 


r < 30 pc 


0.55 ±0.00 


r < 9 pc 


0.48 ±0.01 


r < 3 pc 


0.54 ±0.01 



The bottom panel shows the comparison of the cumu- 
lative mass profiles for the observations and the ensemble 
en3875W12 (A^o = 3875, Wq = 12). This model has been 
chosen because the agreement with the observations is very 
good in the core (r < 3 pc) and with respect to the final 
total mass contained within a radius of r = 30 pc from the 
cluster center. In the range 3 pc < r < 18 pc pc the fit is 
fairly bad, i.e. not consistent within 1(t„i. A better fit of the 
core region cannot be obtained for Wq — 12. For Wq — 3 
the models dissolve so fast that a good fit with respect to 
the inner cumulative mass profile or the final total mass 
has not been obtained at all. In Appendix |B] we show the 
best fits for the present-day total mass within r = 30 pc 
from the center, which, however, are not agreeing well with 
respect to the inner cumulative mass profile. 
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Fig. 4. Comparison of the present-day cumulative mass 
profiles of observations and models. The thick solid (red) 
lines in the middle for the models represent the mean which 
is ensemble-averaged over 15 runs. The standard error is 
shown as a filled area. The dashed (blue) line is calculated 
from the observational data. Top: The ensemble en3000W6. 
Middle: The ensemble en3000W9. Bottom: The ensemble 
en3875W12. 



5.2. Past evolution and present-day state 

The best fitting ensemble of models (en3000W6) with re- 
spect to the inner cumulative mass profile is an ensemble 
average over 15 runs with different random number seeds. 
The initial and final parameters of the ensemble en3000W6 
are given in Table[3j The final parameters are obtained after 
625 Myrs of evolution of the initial model and correspond 
to the present-day state. 

In the following discussion, we look at stars within dis- 
tances of r = 30 pc from the cluster center, within r = 9 
pc (Jacobi radius) and at stars within r = 3 pc (core). 

Figure [5] shows the evolution of the particle number 
within r = 30 pc from the cluster center. It can be seen 
how the 15 cluster models of the ensemble en3000W6 dis- 
solve as time proceeds. The mean final particle number and 
the standard deviation for a single run is shown as a (red) 
dot with errorbars at t = 625 Myrs. The standard devia- 
tion is rather high due the non-Markovian evolution of the 
models which produces rms scatter between the individual 
models of an ensemble. The average particle loss rate until 
t = 625 Myr is dN/dt = — 3.6/Myr. The average mass loss 




Fig. 5. Particle number within r = 30 pc from the cluster 
center as a function of time for all 15 runs of the ensemble 
en3000W6. One can see the non-Markovian process of star 
cluster dissolution. The mean final particle number and the 
standard deviation for a single run is shown as a (red) dot 
with vertical errorbars sX t = 625 Myrs. The future evolu- 
tion of the mean particle number within r = 30 pc is shown 
in the inset. There the times when the cluster has reached 
20, 10 and 5% of its initial particle number are shown as 
(red) dots. 
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Fig. 6. The present-day state of one run of the ensemble 
of models en3000W6 after 625 Myrs of evolution. The tidal 
tails have reached a length of 800 pc. 



rate until t = 625 Myr is dM/dt = -2.2 Mg/Myr. The 
inset of Figure [5] is explained in Section |5.3| 

Figure |6] shows the present-day state of the cluster with 
its tidal tails for one run of the ensemble en3000W6. The 
tidal tails have reached a length of 800 pc after 625 Myrs 
of evolution. However, the kicked out NSs and BHs (and a 
few WDs) have hurried ahead along the cluster orbit during 
the evolution and reside dX t = 625 Myrs spread out along 
the orbit well beyond the tips of the tidal tails at 800 pc. 

Table |4] shows the values of the 3-dimensional veloc- 
ity dispersion ctsd for the ensemble en3000W6 within 30, 
9 and 3 pc from the cluster center. The three values are 
consistent with the fact that the velocity dispersion in the 
core is higher than in the halo within the Jacobi radius and 
rises again outside of the Jacobi radius. The values can be 
compared with those given in Roser et al. (2011). 

5.2.1. Stellar mass functions 

Figure I?] shows the average reahzation of the Kroupa (2001) 
IMF of the ensemble en3000W6 with lowest mass mi — 
0.08 A/0 and initial particle number TVq = 3000. Shown is 
the mean value of \og{N) from the ensemble averaging over 
15 runs. The errors are 2am standard errors. The analytical 
slopes and the transition at 0.5 Mq are also shown as lines. 
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Fig. 7. The average realization of the Kroupa (2001) IMF 
of the ensemble en3000W6 with lowest mass to; — 0.08 Mq, 
highest mass nih = 100.00 Mq and initial particle number 
A'o = 3000. Shown is the mean value of log(7V) from the en- 
semble averaging over 15 runs. The Icjm standard errors are 
shown. The analytical slopes and the transition at 0.5 Mq 
are also shown as lines. Note that the depletion in the lowest 
mass bin is due to the binning. 



Figure [8] shows a comparison of the PDMFs between ob- 
servations and ensemble en3000W6. The top panel shows 
the stars within a radius of r = 30 pc from the center; the 
middle panel shows the stars within a radius of r = 9 pc 
from the center, i.e. within the Jacobi radius; the bottom 
panel shows the stars within a radius of r = 3 pc from 
the center (core). Shown is the mean value of log(iV) from 
the ensemble averaging over 15 runs. The errors are Idm 
standard errors. The blue bars are derived from the obser- 
vational data with iV — 724 stars within 30 pc and lowest 
mass rai = 0.116 Mq while the red bars are derived from 
the ensemble of models with (iV) — 736 stars within 30 
pc and lowest mass mi = 0.08 Mq. Note that these two 
particle numbers N of model and observations are indeed 
realized in the top panel. This cannot easily be seen due 
to an "optical illusion" : The deviations between ensemble 
and observation carry more weight at high values of log(iV). 
Furthermore, note that the completeness limit of the obser- 
vations is around 0.25 M0 (Roser et al. 2011). 

The following details can be seen in Figure |8] In all 
panels the upper main sequence (m > 1 Mq) is overabun- 
dant in the observations. The reason may be simply that 
the IMF of our models is wrong . Also , unresolved binaries 
could be the reason (see Section 5.2.4). Binaries are not re- 
solved below 3", i.e. 133 astronomical units at a distance 
of 46 pc (Roser, priv. comm.). In the upper two panels the 



PDMF of the ensemble is nearly flat at the low mass end 
while the observed PDMF decreases below the complete- 
ness limit of the observations. In the upper two panels, a 
dip in the PDMFs between m = 0.5 — 0.6 Mq can be seen. 
This dip in the model coincides with the location of the 
Wielen dip (Kroupa et al. 1990) which is related to the 
shape of the mass- luminosity relation (Kroupa 2002). The 
stellar evolution in our models is different for m < 0.7 Mq 
and TO > 0.7 Mq (Hurley et al. 2000). 

Figure 9] shows the cumulative stellar mass function for 
the ensemble en3000W6 for stars within a radius of r = 9 
pc from the cluster center, i.e. within the Jacobi radius. 
Shown is the sum of all masses which are larger than a 
given mass. The histogram includes main sequence stars, 
giants and WDs. Stellar mass BHs are excluded and NSs do 
not occur. Up to the completeness limit of the observations 
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Fig. 8. Comparison of the PDMFs between observations 
and ensemble of models en3000W6. Top panel: Stars within 
a radius of 30 pc from the center; Middle panel: Stars within 
a radius of 9 pc from the center (Jacobi radius); Bottom 
panel: Stars within a radius of 3 pc from the center (core) . 
Shown is the mean value of log(iV) from the ensemble av- 
eraging over 15 runs. The errors are \am standard errors. 
Note that the blue bars are derived from the observational 
data with N = 724 stars and lowest mass to/ = 0.116 Mq 
while the red bars are derived from the ensemble of models 
with (iV/) = 736 stars and lowest mass mi = 0.08 Mq. The 
onset of incompleteness of the observations is shown as a 
vertical line in the top panel. 



there is a discrepancy of 10 % which leaves room for further 
optimization in future models. 

5.2.2. Stellar evolution 

Table [5] shows the number counts of stellar types for the 
best fitting ensemble of models (en3000W6) for all 3000 
stars (initial and final), for stars within a radius of 30 pc 
(final), 9 pc (final inside Jacobi radius) and 3 pc (final inside 
core). The stellar types are the same as in the classification 
in the beginning of Section 4 of Hurley et al. (2000). The 
standard errors from the ensemble-averaging over 15 runs 
are also given. We find 1 ± stellar mass BHs within r — 
3 pc. Almost all stellar mass BHs and all NSs have been 
kicked out due to supernova kicks and are spread out along 
the cluster orbit. We find 27 ± 2, 22 ± 2 and 9 ± 1 carbon- 
oxygen WDs within r = 30, 9 and 3 pc, respectively. These 
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Fig. 9. Present-day cumulative mass function for the en- 
semble en3000W6 for r < 9 pc from the cluster center. 
Shown is the sum of all masses which are larger than a 
given mass. The errors are Idm standard errors. The com- 
pleteness limit of the observations is shown as a vertical 
solid line. 
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Fig. 10. Present-day synthetic Hertzsprung-Russell dia- 
gram for one run of the evolved model en3000W6. Shown 
are stars within a radius of r = 30 pc from the cluster cen- 
ter. The final particle number within r = 30 pc for this run 
is Nf ^ 836. 



numbers can be compared to the observations (Schilbach & 
Roser 2011). 

Figure [lO] shows a present-day synthetic Hertzsprung- 
Russell diagram for one run of the evolved ensemble 
en3000W6. Shown are stars within a radius of r = 30 pc 
from the cluster center. Note that the synthetic diagram in 
Madsen (2003) is more suited for a direct comparison with 
observations since he transforms effective temperatures to 
B — V colors using color tables and bolometric magnitudes 
to V magnitudes. In addition, he introduces a small scatter 
in both V and B — V. We cannot directly compare to the 
observations since the data in B — ^ are missing for the 
data set used in the present study. 

Figure pT] show the comparison of the present-day lumi- 
nosity functions (PDLFs) of the ensemble en3000W6 and 
the observations by Roser et al. 2011. In the top, middle 
and bottom panels stars within a radius of 30, 9 (Jacobi 



Table 5. Number counts of stellar types for the best fitting 
ensemble of models (en3000W6) for all 3000 stars (initial 
and final), for stars within a radius of 30 pc (final), 9 pc 
(final inside Jacobi radius) and 3 pc (final inside core). The 
stellar number counts are averaged over the 15 runs of the 
ensemble. The stellar types are the same as in the classifi- 
cation in the beginning of Section 4 of Hurley et al. (2000). 



# Type Name N ± a,,, 

1 Total (initial, all stars) 

Low main sequence m < 0.7 Mq 2567 ± 5 

1 Main sequence m > 0.7 M© 432 ± 5 

2 Total (final, all stars) 

Low main sequence m < 0.7 Mq 2567 ± 5 

1 Main sequence m > 0.7 Mq 348 ± 5 
4 Core helium burning 7 ± 

11 Carbon-oxygen WDs 48 ± 2 

12 Oxygen-Neon WDs 7 ± 1 

13 NSs 13 ± 1 

14 BHs 5 ± 

3 r < 30 pc (final) 

Low main sequence m < 0.7 Mq 565 ± 39 

1 Main sequence m > 0.7 Mq 138 ± 10 
4 Core helium burning 4 ± 
11 Carbon-oxygen WDs 27 ± 2 
14 BHs 1±0 

4 r < 9 pc (final, inside Jacobi radius) 

Low main sequence m < 0.7 Mq 343 ± 32 

1 Main sequence m > 0.7 Mq 107 ± 10 
4 Core helium burning 4 ± 
11 Carbon-oxygen WDs 22 ± 2 
14 BHs 1±0 

5 r < 3 pc (final, inside core) 

Low main sequence m < 0.7 Mq 70 ± 8 

1 Main sequence m > 0.7 Mq 42 ± 5 
4 Core helium burning 2 ± 
11 Carbon-oxygen WDs 9 ± 1 
14 BHs 1 ± 



radius) and 3 pc (core) from the cluster center have been 
used in the statistics. Shown is the mean value of log(7V) 
as a function of the 2MASS Kg band magnitude (observa- 
tions) and K band magnitude (model) from the ensemble 
averaging over 15 runs. The errors are lam standard errors. 

In order to obtain the absolute K band magnitudes 
from the theoretical luminosities given by nbody6 we ap- 
plied the compilation of infrared bolometric corrections 
given in Appendix A of Muench (2002) with the calibra- 
tion mboi,0 = 4.75. The applied conversion formula is 

Mk = -2.51ogio(L/L0) + Mboi,0 - BCk{T,s), (13) 

where Mk, L/Lq, Afboi.o and BCk are the absolute K 
band magnitudes, the luminosity in units of the solar lu- 
minosity, the bolometric magnitude of the Sun and the K 
band bolometric correction, respectively. 

In Figure [TT] there are giants in the bin centered at 
Mk = —1.75 of the observational LF. For these luminous 
objects it is hard to obtain the correct observational Mk 
due to overexposure (S. Roeser and E. Schilbach, priv. 
comm.). On the other hand, in the theoretical LF the 
WDs with logiQ(rcff) = 4.0 - 4.4 populate the faint end 
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Fig. 11. PDLFs in the K band for the ensemble en3000W6 
and the observations (Roser et al. 2011). For explanations 
see the text. Shown is the mean/measured value of log(A^) 
as a function of the absolute Kg (observations) / K (model) 
band magnitude. The errors are \<Jm standard errors from 
the ensemble averaging. 



{Mk > 10). Between 1 < Mk < 5 it can be seen that the 
observational Mj^^s are overabundant as compared to the 
model Mk's. 



5.2.3. Mass segregation 

A first unmistakable sign of mass segregation is the in- 
creased mean mass within the core in Table |3] as compared 
to the mean mass within the Jacobi radius. 

Figure [12] shows the average cumulative mass for obser- 
vation and ensemble en3000W6. The average cumulative 
mass is calculated only for stars with mass m > 0.116 Mq. 
Also shown is the dotted (magenta) line for the ensemble 
enb2125PL with 33% primordial binaries (see Section 5.2.4 
below) . The Figure clearly shows the segregation of masses 
within the Jacobi radius. 



Figure [T3| shows the result of a detailed analysis of mass 
segregation based on the minimum spanning tree (MST) 
method A^^gx developed by Olczak, Spurzem & Henning 
(2011). They define a measure of mass segregation Fmst, 
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Fig. 12. Average cumulative mass for observations and en- 
sembles en3000W6 and enb2125PL. The average cumu- 
lative mass is calculated only for stars with mass m > 

0.116 Mq. 
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where Ci are the lengths of the n MST edges. The super- 
script "ref" refers to a sample of n -I- 1 random stars from 
the entire population while the superscript "mass" refers 
to an equal-sized sample of the most massive stars. Note 
that 7mst has the dimension of length while Fmst is a di- 
mensionless measure. Fmst • (AFmst) > 1 implies mass 
segregation of the n+1 most massive stars with significance 
k<7. 

The top panel of Figure [13] shows the analysis of the 
evolved (for 625 Myrs) ensemble en3000W6 while the bot- 
tom panel shows an analysis of the observational data 
of Roser et al. (2011). Each single plot contains on the 
left hand-side the "cumulative" Fmst (points) for the 5 
(red), 10 (green), 20 (blue), 50 (magenta), 100 (cyan), 200 
(brown), 500 (orange), and all most massive stars (black). 
On the right-hand side the "differential" Fmst (points) is 
shown for the 5 (red), 6-10 (green), 11-20 (blue), 21-50 (ma- 
genta), 51-100 (cyan), 101-200 (brown), 201-500 (orange), 
and 501-all most massive stars (black). The error bars mark 
the la, 2(7, and 3(J uncertainties. A value of one (dashed 
line) marks the unsegregated state. The higher the values, 
the more segregated the given mass group is (compared to 
a set of random groups) . The results of this analysis are the 
following: 

1. The ensemble average of the numerical models resem- 
bles the observational data very well: besides a larger 
deviation of the 4th bin the agreement is excellent. 

2. Both observations and simulations show a significant 
degree of mass segregation, mostly above the 3cr level 
for the cumulative analysis. 

3. There is a clear signature of "inverse mass segregation" : 
Tmst < 1 for the last bin in the differential plot (i.e. 
the ~ 200 least massive stars). This demonstrates the 
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Fig. 13. Comparison of present-day mass segregation be- 
tween models and observations. Top figure: Analysis of the 
evolved (for 625 Myrs) ensemble en3000W6. Bottom figure: 
Analysis of the observational data of Roser et al. (2011). 
Each single plot contains on the left hand-side the "cu- 
mulative" Tmst (points) for the 5 (red), 10 (green), 20 
(blue), 50 (magenta), 100 (cyan), 200 (brown), 500 (or- 
ange), and all most massive stars (black). The right-hand 
side shows the "differential" T mst (points) for the 5 (red) , 
6-10 (green), 11-20 (blue), 21-50 (magenta), 51-100 (cyan), 
101-200 (brown), 201-500 (orange), and 501-all most mas- 
sive stars (black). The error bars mark the 1, 2, and Scr 
uncertainties. A value of one marks the unsegregated state. 
The higher the values, the more segregated the given mass 
group is (compared to a set of random groups) . 



removal of the lowest mass members from the inner clus- 
ter parts due to dynamical mass segregation. 
The signature of mass segregation agrees fairly well with 
a moderately [S = 0.3: cf. Subr, Kroupa & Baumgardt 
2008) mass segregated model star cluster of 1000 parti- 
cles as shown in the middle panel of Fig. 4 in Olczak, 
Spurzem & Henning (2011). 




Fig. 14. Comparison of the present-day cumulative mass 
profiles of observations and models with primordial bina- 
ries. The thick solid (red) lines is derived from the en- 
semble en2250PL without binaries, the dashed (green) line 
is derived from the ensemble enb2250PL with 33% pri- 
mordial binaries and the dotted (magenta) line is derived 
from the best fitting ensemble with 33% primordial binaries 
enb2125PL. The dashed (blue) line represents the observa- 
tions. The standard errors are shown as filled areas. 

Table 6. Parameters of the ensemble enb2125PL with 33% 
primordial binaries. Notes: (1) The final parameters are ob- 
tained after 625 Myrs of evolution, i.e. at the present time. 
(2) The final "observed" particle number is the number of 
stars with m > 0.116 Mq. (3) The highest mass is obtained 
by excluding BHs. 



# 


Quantity 




Value 


1 


Initial model (initial) 








Particle number A'^o 




2125 




Total mass (Mo) 




1230 Mq 




Jacobi radius rj 




14.5 pc 


2 


r < 30 pc (final) 








Particle number (A^f) ± o"™ 




619 ±58 




"Observed" part. numb. (At)obs 


± Om 


542 ± 51 




Lowest mass m; 




0.08 Mq 




Highest mass nih ± Om 




2.86 ± 0.02 Mq 




Mean mass (m) 




0.606 Mq 




Total mass Mf ± am 




375 ± 35 Mq 


3 


r < 9 pc (final, inside Jacobi radius) 






Particle number (Af) ± (Jm 




403 ± 49 




"Observed" part. numb. (At)obs 


± Om 


358 ± 44 




Mean mass (m) 




0.679 Mq 




Total mass Mf ± am 




274 ± 33 Mq 


4 


r < 3 pc (final, inside core) 








Particle number (Af) ± am 




127 ± 21 




"Observed" part. numb. (Af)obs 


± O-m 


116 ± 20 




Mean mass (m) 




0.882 Mq 




Total mass M{ ± am 




112 ± 19 Mq 



5.2.4. Primordial binaries 

The recent work of SoUima et al. (2009) showed that open 
clusters may contain between « 35 — 70% binaries. We 
have run the best fitting ensemble en2250PL of Plummer 
models again with a prim ordial binary fraction = 33% 



{fsb = 50% from Eq. |10[ ) in order to investigate the influ- 
ence of primordial binaries. The ensemble has been called 



enb2250PL. The third letter "b" denotes primordial bina- 
ries and the following number is the number of stars (i.e. 
not the number of systems). Since this ensemble does not 
fit well to the inner observational cumulative mass profile 
we further reduced the initial particle number and finally 
found a best fitting model enb2125PL. 
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Fig. 15. PDLFs in the K band for the ensemble enb2125PL 
where all binaries are either resolved or unresolved. Top 
panel: Stars within a radius of 30 pc from the center; Middle 
panel: Stars within a radius of 9 pc from the center (Jacobi 
radius); Bottom panel: Stars within a radius of 3 pc from 
the center (core). Shown is the mean value of log(7V) as a 
function of the absolute K band magnitude. 



Fig. 16. Comparison of the PDMFs between ensembles 
enb2250PL, en2250PL and enb2125PL. Middle panel: Stars 
within a radius of 9 pc from the center (Jacobi radius); 
Bottom panel: Stars within a radius of 3 pc from the center 
(core). Shown is the mean value of log(A^) from the ensem- 
ble averaging over 15 runs. The errors are l<7m standard 
errors. 



Figure [14] shows the comparison of the present-day cu- 
mulative mass profiles of observations and models with 
primordial binaries. The thick solid (red) lines is derived 
from the ensemble en2250PL without binaries, the dashed 
(green) line is derived from the ensemble enb2250PL with 
33% primordial binaries and the dotted (magenta) line is 
derived from the best fitting ensemble with 33% primor- 
dial binaries enb2125PL. The dashed (blue) line repre- 
sents the observations. The standard errors are shown as 
filled areas. It can be seen that the profile of the ensemble 
enb2250PL is steeper than that of the ensemble en2250PL 
at small radii. The ratio of half-mass to Jacobi radius for 
the model enb2125PL is given by rh/rj = 0.18 ± 0.01. 
The field star contamination for the model enb2125PL is 
|M^odci-Mobs|/Mobs = 12% (3 pc), 0% (9 pc), 13% (18 pc), 
20% (30 pc). 

Figure [15] shows the PDLFs in the K band for the en- 
semble enb2125PL where all binaries are either resolved or 
unresolved. For the unresolved binaries we transformed for 
each star the bolometric luminosities given by nbody6 to 
absolute K band magnitudes using the bolometric correc- 
tions of Muench (2002) and then back to K band lumi- 
nosities using the absolute K-band magnitude of the Sun 



Mk,o = 3.33. Then we added the K band luminosities of 
both binary components and transformed back to absolute 
K band magnitudes which are shown in Figure |15[ With 
the present parameters of the distribution of binaries unre- 
solved binaries produce a difference at the faint end of the 
LF. However, our statistics of the binary parameters may 
be not fully realistic. 

Figure [16] shows a comparison of the PDMFs between 
ensembles enb2250PL, en2250PL and enb2125PL. It can be 
seen that the ensembles en2250PL and enb2125PL are de- 
pleted of high-mass stars as compared with the ensemble 
enb2250PL for r < 30 pc and r < 9 pc. Furthermore, the 
comparison of ensembles en2250PL and enb2125PL shows 
that the ensemble enb2125PL is depleted of low-mass stars 
as compared with the ensemble en2250PL while high-mass 
stars are overabundant in the ensemble enb2125PL as com- 
pared with the ensemble en2250PL. 

The initial and final parameters of the ensemble 
enb2125PL are given in Table ]6] The final parameters are 
obtained after 625 Myrs of evolution of the initial model 
and correspond to the present-day state. 

Table [7] shows the number counts of stellar types that 
are the component of a binary for the ensemble enb2125PL 
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Table 7. Number counts of stellar types that are the com- 
ponent of a binary for the ensemble enb2125PL for all 2125 
stars (final), for stars within a radius of 30 pc (final), 9 
pc (final inside Jacobi radius) and 3 pc (final inside core). 
The stellar number counts are averaged over the 15 runs 
of the ensemble. The stellar types are the same as in the 
classification in the beginning of Section 4 of Hurley et al. 
(2000). 



# 


Type 


Name 


N ±o^ 


1 


Total 


(final, all stars) 









Low main sequence m < 0.7 M© 


964 ±4 




1 


Main sequence m > 0.7 Mq 


18 ±2 




11 


Carbon-oxygen WDs 


12 ± 1 


2 


r < 30 pc 


(final) 









Low main sequence m < 0.7 Mq 


513 ±50 




1 


Main sequence m > 0.7 Mq 


9± 1 


3 


r < 9 pc 


(final, inside Jacobi radius) 









Low main sequence m < 0.7 Mq 


346 ± 44 




1 


Main sequence m > 0.7 Mq 


8± 1 


4 


r < 3 pc 


(final, inside core) 









Low main sequence m < 0.7 Mq 


111 ± 19 




1 


Main sequence m > 0.7 Mq 


7± 1 



for all 2125 stars (final), for stars within a radius of 30 pc 
(final), 9 pc (final inside Jacobi radius) and 3 pc (final inside 
core). The stellar number counts are averaged over the 15 
runs of the ensemble. The stellar types are the same as in 
the classification in the beginning of Section 4 of Hurley 
et al. (2000). There are no BHs or NSs in binaries. Also, 
there are no WDs in binaries within a radius of 30 pc from 
the cluster center. They are all kicked out due to binary 
kicks. Since WDs in binaries are observed in the Hyades 
(Schilbach & Roser 2011) they may not be subject to strong 
kicks. 



5.3. Future evolution 

The future evolution of the ensemble en3000W6 is shown 
in the inset of Figure [5 There the times when the cluster 
has reached 20, 10 and 5% of its initial particle number 
are shown (magenta) with horizontal error bars. The times 
are <20 = 695 Myrs, tio = 875 Myrs and is = 1029 Myrs. 
The latter time may be taken as an approximate dissolution 
time. 



6. Discussion 

The present-day properties of the Hyades can be well repro- 
duced by a standard King or Plummer initial model when 
choosing appropriate initial conditions. We have found a 
model (en3000W6) which imitates the cumulative mass pro- 
file within the Jacobi radius of the observed Hyades spatial 
distribution very well. Also, the derived PDMF and the K 
band LP show a good agreement between observations and 
model. 

From our models we can make the following statements 
about the Hyades: 

— The tidal tails of the Hyades must have a length of w 800 
pc today. 



— Roche-lobe filling Wq = 6 — 9 single-star King initial 
models (King 1966) with the initial particle number 
A'o — 3000 provide a good fit to the inner present-day 
cumulative mass profile. 

— The best-fit single-star roche-lobe filling King model has 
an initial mass of Mq = 1721 Mq, an initial Jacobi 
radius of rj = 16.2 pc and an average mass loss rate of 
dM/dt = 2.2 Mo/Myr. 

— A reasonably good fit is obtained with a Plummer model 
with 33% primordial binaries and a ratio of half-mass to 
Jacobi radius of Vh/rj — 0.18. Here the initial particle 
number and mass are Nq — 2125 and Mq — 1230 Af©, 
respectively. The average mass loss rate is dM jdt — 
1.4 M0/Myr. The observed average cumulative mass is 
reproduced very well. 

— Mass segregation is clearly detected in both observa- 
tions and models of the Hyades cluster. 

— The Hyades contain only 1 ± stellar mass black holes. 
The probability is very high that nearly all NSs and 
BHs are kicked out due to supernova kicks. 

— The number of WDs critically depends on the kick ve- 
locity which is adopted for WD kicks. This kick veloc- 
ity is not yet well constrained. Assuming a mean kick 
velocity of « 8 km/s (which is comparable with the es- 
cape velocity from the cluster center) we obtain 27 ± 2 
carbon-oxygen WDs within a radius r < 30 pc from the 
cluster center, 22 ± 2 carbon-oxygen WDs within the 
Jacobi radius (r < 9 pc) and 9 ± 1 carbon-oxygen white 
dwarfs within the core (r < 3 pc). For the comparison 
to the observations we refer to the article by Schilbach 
& Roser (2011). We do not find WDs in binaries within 
r < 30 pc in our model. They are all kicked out in the 
form of binary kicks. However, some WDs in binaries are 
observed in the Hyades (Schilbach & Roser 2011). This 
suggests that WDs in binary systems are not subject to 
strong kicks. 

— The inclusion of 33% primordial binaries reduces the 
initial number of Hyades stars A'o by w 5% as compared 
to the model without primordial binaries. 

— Under the asumption that the single-star ensemble 
en3000W6 is a good approximation to reality we find 
that the Hyades will be dissolved up to 5% of its initial 
number of stars in k, 400 Myrs from now. 

The degeneracy of good-fitting models can be quite high 
due to the large dimension of the parameter space, which is 
spanned by the initial particle number A'o , the King param- 
eter Wq ^-nd two extra dimensions: The primordial binary 
fraction /b and the Roche- lobe fiUing factor rgQ%/rj. In this 
work we did not explore in detail the latter two degrees of 
freedom. More simulations with different Roche-lobe fill- 
ing factors and primordial binary fractions are required to 
explore this degeneracy in more detail. 
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Appendix A: nbodyGtidgpu 

The program nbodyGtidgpu is a new direct A^-body code 
for the integration of the A^-body problem in an ana- 
lytic background potential using Graphics Processing Units 
(GPUs). The program is essentially based on nbodyG by 
Aarseth (Aarseth 1999, 2003). nbodyGtidgpu is not paral- 
lelized with MPI routines but runs on multi-core processors 
(via OpenMP) with one or two GPU(s). For this purpose, 
special routines in CUDA have been added by Nitadori in 
collaboration with Aarseth. The variant of NBODy6 with 
GPU support has been termed nbody6gpu (Nitadori & 
Aarseth 2011). 

The implementation of different galactic tidal fields into 
NBODy6gpu is the contribution of some of the present au- 
thors and the special feature of nbody6tidgpu. The imple- 
mentation is based on the implementation in the (parallel) 
code nbodyBgc which is in detail described in Ernst, Just 
& Spurzem (2009) and Ernst (2009). Several analytic mod- 
els for the background potential are part of the code, for 
example 

— the three-component Plummer-Kuzmin model of the 
Milky Way described in Section [2] (Miyamoto & Nagai 
1975) 

— power-law models with and without central black hole 

— Dehnen models (Dehnen 1993, Tremaine et al. 1994) 

— Plummer model (Plummer 1911) 

In addition to the equations of motion of the A^-body 
problem, the program solves the equations of motion of 
the star cluster orbit around the Galactic center with an 
eighth-order composition scheme (Yoshida 1990; for the co- 
efficients see McLachlan 1995). The tidal force of the back- 
ground potential acts on all particles in the A^-body system 
and is added as a perturbation to the KS regularization 
(Kustaanheime & Stiefel 1965) of nbody6. 

The program is able to handle the integration of star 
cluster orbits with allmost all eccentricities. For nearly ra- 
dial orbits a time transformation can be switched on in 
order to guarantee an exact orbit integration during the 
pericenter passages. The leapfrog-like symplectic integra- 
tion schemes allow for adaptive time steps if one applies 
a Sundman transformation to the time (Sundman 1912, 
Mikkola & Tanikawa 1999, Preto & Tremaine 1999, Mikkola 
& Aarseth 2002). 

It is also possible to switch on dynamical friction 
with different x functions and different realizations of the 
Coulomb logarithm, i.e. fixed or variable according to Just 
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(A^o = 3375,1^0 = 9). It can be seen that the agreement 
between observations and ensemble of models are rather 
bad with respect to the inner cumulative mass profile, i.e. 
not even consistent within 2am for the largest part of the 
radial range. 



Fig. B.l. Comparison of the cumulative mass profiles of 
observations and ensemble of models. The thick solid (red) 
lines in the middle for the model represent the mean which 
is ensemble-averaged over 15 runs. The standard error of 
the mean for the ensemble is shown as a filled area. The 
dashed (blue) line is calculated from the observational data. 
Top: The ensemble of models en3125W6. Bottom: The en- 
semble of models en3375W9. 



& Penarrubia (2005). Since the symplectic composition 
schemes are by construction suited for Hamiltonian sys- 
tems, the dissipative dynamical friction force requires spe- 
cial attention. It is implemented in nbody6tidgpu with 
an iterative implicit midpoint method (see, e.g., Mikkola & 
Aarseth 2002). 

For the current study, the three- component Plummer- 
Kuzmin model described in Section [2] is used as a back- 
ground potential. The acceleration for the Plummer- 
Kuzmin model is given by the expression 



A 



^ ] ^ I y 



(A.l) 

For the tidal correction to the equations of motion of the 
A'^-body problem, the time derivative of acceleration (jerk) 
of the analytic background potential is needed in order to 
guarantee energy conservation and exactness of the orbit 
integration on a high level. Note that the jerk is given by 
lengthy expressions which depend also on the velocity com- 
ponents Vx, Vy and Vz and which we do not reproduce here. 
However, they are as well implemented in nbodyBtidgpu. 



Appendix B: Bad fits 



Figure |B.1| shows the two best fits for the total mass of 
Hyades members within 30 pc from the cluster center from 
the models within the second parameter grid of Section 4j 
i.e. the grid with AiV = 125. The top panel shows the 
ensemble of models en3125W6 (Nq = 3125, Wq = 6) while 
the bottom panel shows the ensemble of models en3375W9 
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